%**************************************************************************
% Business-as-Usual and Orbital Debris Path
% Code for the Figures
% A. Bongers and J. L. Torres (2025)
% **************************************************************************
clear all
close all

nu1=33.33333;

% Clean space
file0='soldise00.xlsx';
y0=xlsread(file0,'y','c4:c204');
i0=xlsread(file0,"ik","c4:c204");
h0=xlsread(file0,"is","c4:c204");
k0=xlsread(file0,"k","c4:c204");
s0=xlsread(file0,"s","c4:c204");
S0=xlsread(file0,"ss","c4:c204");
H0=xlsread(file0,"hh","c4:c204");
L0=xlsread(file0,"ll","c4:c204");
c0=xlsread(file0,"c","c4:c204");

% Central planner
file1='soldise01.xlsx';
y1=xlsread(file1,'y','c4:c204');
i1=xlsread(file1,"ik","c4:c204");
h1=xlsread(file1,"is","c4:c204");
k1=xlsread(file1,"k","c4:c204");
s1=xlsread(file1,"s","c4:c204");
S1=xlsread(file1,"ss","c4:c204");
%H1=xlsread(file1,"hh","c4:c204");
L1=xlsread(file1,"ll","c4:c204");
D1=xlsread(file1,"dd","c4:c204");
X1=xlsread(file1,"xx","c4:c204");
Z1=xlsread(file1,"z","c4:c204");
W1=xlsread(file1,"w","c4:c204");
F1=xlsread(file1,"ff","c4:c204");
eds1=xlsread(file1,"dso","c4:c204");
erb1=xlsread(file1,"rb","c4:c204");
ecol1=xlsread(file1,"col","c4:c204");
edsb1=xlsread(file1,"dsb","c4:c204");
erbb1=xlsread(file1,"rbb","c4:c204");
efrag1=xlsread(file1,"efrag","c4:c204");
ecolds1=xlsread(file1,"dsd","c4:c204");
ecolrb1=xlsread(file1,"rbd","c4:c204");
c1=xlsread(file1,"c","c4:c204");

% Decentralized economy
file2='decentralized.xlsx';
y2=xlsread(file2,'y','c4:c204');
i2=xlsread(file2,"ik","c4:c204");
h2=xlsread(file2,"is","c4:c204");
k2=xlsread(file2,"k","c4:c204");
s2=xlsread(file2,"s","c4:c204");
S2=xlsread(file2,"ss","b4:b204");
%H2=xlsread(file2,"hh","c4:c204");
L2=xlsread(file2,"ll","b4:b204");
D2=xlsread(file2,"dd","b4:b204");
X2=xlsread(file2,"xx","b4:b204");
Z2=xlsread(file2,"z","b4:b204");
W2=xlsread(file2,"w","b4:b204");
F2=xlsread(file2,"ff","b4:b204");
eds2=xlsread(file2,"dso_temp","b4:b204");
erb2=xlsread(file2,"rb_temp","b4:b204");
ecol2=xlsread(file2,"ex_temp","b4:b204");
edsb2=xlsread(file2,"edsb_temp","b4:b204");
erbb2=xlsread(file2,"erbb_temp","b4:b204");
efrag2=xlsread(file2,"efrag_temp","b4:b204");
ecolds2=xlsread(file2,"eds_temp","b4:b204");
ecolrb2=xlsread(file2,"erb_temp","b4:b204");
c2=xlsread(file2,"c","c4:c204");

T=125;
Periods=2023:2023+T-1;

% Figure 1
figure
subplot(2,2,1)
plot(Periods,D1(1:T),Periods,D2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of debris')
title('(a) Number of orbital debris > 1 cm')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,2)
plot(Periods,F1(1:T),Periods,F2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(b) Number of fragments > 1 cm')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,3)
plot(Periods,W1(1:T),Periods,W2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of satellites')
title('(c) Derelict satellites')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,4)
plot(Periods,Z1(1:T),Periods,Z2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Rocket bodies')
title('(d) Rocket bodies')
legend('Central planner','Laissez-faire')
grid on

% Figure 2
figure
subplot(2,2,1)
plot(Periods,S1(1:T),Periods,S2(1:T),'--',Periods,S0(1:T),':',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number')
title('(a) Satellites')
legend('Central planner','Laissez-faire','Clean space')
grid on
subplot(2,2,2)
plot(Periods,L1(1:T),Periods,L2(1:T),'--',Periods,L0(1:T),':',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number')
title('(b) Launches')
legend('Central planner','Laissez-faire','Clean space')
grid on
subplot(2,2,3)
plot(Periods,X1(1:T),Periods,X2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number')
title('(c) Satellites destroyed')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,4)
plot(Periods,1.25*10^-10*D1(1:T),Periods,1.25*10^-10*D2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Probability of collision')
title('(d) The Kessler syndrome')
legend('Central planner','Laissez-faire')
grid on

% Figure 3
figure
subplot(2,2,1)
plot(Periods,(1+nu1)*ecol1(1:T),Periods,(1+nu1)*ecol2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(a) Satellite collisions')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,3)
plot(Periods,(1+nu1)*eds1(1:T),Periods,(1+nu1)*eds2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(c) Derelict satellites')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,4)
plot(Periods,(1+nu1)*erb1(1:T),Periods,(1+nu1)*erb2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(d) Rocket bodies')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,2)
plot(Periods,(1+nu1)*efrag1(1:T),Periods,(1+nu1)*efrag2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(b) Launches')
legend('Central planner','Laissez-faire')
grid on

% Figure 4
figure
subplot(2,2,1)
plot(Periods,(1+nu1)*edsb1(1:T),Periods,(1+nu1)*edsb2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(a) Derelict satellite breakups')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,2)
plot(Periods,(1+nu1)*erbb1(1:T),Periods,(1+nu1)*erbb2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(b) Rocket bodies breakups')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,3)
plot(Periods,(1+nu1)*ecolds1(1:T),Periods,(1+nu1)*ecolds2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(c) Derelict satellite collisions')
legend('Central planner','Laissez-faire')
grid on
subplot(2,2,4)
plot(Periods,(1+nu1)*ecolrb1(1:T),Periods,(1+nu1)*ecolrb2(1:T),'--',LineWidth=2) 
xlabel('Periods (years)')
ylabel('Number of fragments')
title('(d) Rocket bodies collisions')
legend('Central planner','Laissez-faire')
grid on
